Before we even start, it is of crucial importance to first understand deeply linear and logistic regression. This is because (as we will see in later discussion) they are the very basic components in a general neural network model.
A linear model can be written in a matrix form
\[ Y = X\beta + \epsilon, \]
where \(Y\) is the output or label vector of length \(N\) (number of observations), \(X\) is the input feature matrix (referred to as the design matrix in statistics) with dimension \(N\) by \(P\) (number of features), \(\beta\) is the model weights/coefficients (a column vector of length \(P\)) for which we’d like to solve, \(\epsilon\) is the model residual or error vector.
The classical way to solve for \(\beta\) is ordinary least squares. The idea of OLS is find out the weights that minimize the mean of squared model errors:
\[ \begin{aligned} \mbox{mse (loss)} &= \frac{1}{N}\sum_i\epsilon^2_i \\ &= \frac{1}{N}\sum_i(y_i - \beta x_i)^2, \end{aligned} \]
where \(y_i\) and \(x_i\) is the \(i\)-th observation. The first-order condition (requiring that the first-order derivative w.r.t. weights are zero) gives us the OLS solution for model weights \(\beta\) analytically (in matrix notation):
\[ \begin{equation} \label{eq:ols} \hat{\beta} = (X'X)^{-1}X'Y. \end{equation} \]
import numpy as np
from numpy.linalg import inv
def ols(X, y):
return inv(X.T.dot(X)).dot(X.T).dot(y)Let’s consider a toy model with only one non-constant feature:
\[ y_i = 6 + 4x_i + \epsilon_i. \]
In this model the outcome \(y\) is determined by a bias term plus a single variable \(x\), with an independently distributed noise term \(\epsilon \sim \mbox{Normal}(0, 1)\).
Create some random data generated from this model:
# Create toy example.
np.random.seed(777)
N = 1000
X = np.stack([np.ones(N), np.random.normal(size=N)], axis=1)
beta = np.array([6, 4], dtype=np.float32) # True model weights.
e = np.random.normal(size=N)
y = X.dot(beta) + e
print(X[:10])from plotly.offline import plot
import plotly.graph_objs as go
def plot_offline(data, layout, ofile):
p = plot({"data": data, "layout": layout},
filename=ofile, auto_open=False, include_plotlyjs=False)
return p
ofile = "plots/toy_reg.html"
p = plot_offline(
ofile=ofile,
data=[go.Scatter(x=X[:,1], y=y, mode="markers")],
layout=go.Layout(title="Data Generated from Toy Model",
xaxis=dict(title="x"),
yaxis=dict(title="y")))Without consideration of the noise term, the OLS estimator will solve precisely for the true model weights:
Of course in the real world the noise term cannot be determined and usually the feature alone cannot explain entirely the variation in the outcome. This means that what we actually estimate will be the expected value of model weights:
We can check the result from sklearn:
from sklearn.linear_model import LinearRegression
# We set fit_intercept=False since our design matrix already contains intercept.
print(LinearRegression(fit_intercept=False).fit(X, y).coef_)By the Law of Large Number and Central Limit Theorem, the OLS estimator will converge in probability to the true model weights and distributed asymptotically Normal in large sample.1
The issue of the above approach is that equation \(\eqref{eq:ols}\) is not numerically stable when it comes to large-scale application where we may have lots of observations and lots of features. One very useful solution to solve the estimator numerically in large-scale application is the gradient descent approach.
Instead of solving the first-order condition analytically, we can do it numerically. Gradient descent is a 1st-order optimization technique to find local optimum of a given function.
In the model training exercise our target function is the loss so the optimization problem is:
\[ \operatorname*{argmin}_\beta \mbox{Loss} = \frac{1}{N}\sum_i(y_i - \beta x_i)^2. \]
That is, we’d like to figure out model weights that minimize the loss which is defined by the mean squared errors when the model is a regression model.
The idea of gradient descent is to
Let’s use the toy example to actually implement a gradent descent optimizer from scratch. First we re-write the loss function explicitly with our setup of one coefficient with a constant (\(\beta = [\beta_0, \beta_1]\)):
\[ \mbox{Loss} = \frac{1}{N}\sum_i\big[y_i-(\beta_0 + \beta_1x_i)\big]^2. \]
Now the gradient (or equivalently the 1st-order derivative) w.r.t. to weights will be:
\[ \begin{aligned} \frac{\partial\mbox{Loss}}{\partial\beta_0} &= - \frac{2}{N}\sum_i \big[ y_i - (\beta_0 + \beta_1x_i) \big], \\ \frac{\partial\mbox{Loss}}{\partial\beta_1} &= - \frac{2}{N}\sum_i \big[ y_i - (\beta_0 + \beta_1x_i) \big]x_i. \end{aligned} \]
The corresponding python function can be coded as:
If we set the above equations to zero we can solve for \(\beta_0\) and \(\beta_1\) analytically and the solution will be exactly just equation \(\eqref{eq:ols}\). But as we just pointed out it suffers from numerical stability issue.
The minimum implementation of our gradient descent optimizer is just a few lines:
Training Steps 10 | Estimate: [1.17474197 2.12281056]
Training Steps 50 | Estimate: [3.87639887 3.20410217]
Training Steps 100 | Estimate: [5.24246911 3.71384968]
Training Steps 500 | Estimate: [6.01109585 3.96245836]
Training Steps 1000 | Estimate: [6.01132195 3.96247176]
Training Steps 3000 | Estimate: [6.01132196 3.96247175]As we can see the result is approaching to our analytical solution as the number of steps grows.
Learning rate is a hyper-parameter for gradient descent optimizer. The gradient update to our model weights is scaled down by the learning rate to make sure convergence. Too large the learning rate will explode the gradient. Too small the learning rate will slow down the convergence and sometimes result in the optimizer trapped at local sub-optimum.
Let’s re-write our gradient descent optimizer to also track the loss from each training step. And we use the same initialization for a fair comparison.
def loss(X, y, beta):
return ((X.dot(beta) - y)**2).mean()
def gd_optimize(X, y, grad_func, lr=.01, n_step=100):
b = np.array([0.0, 0.0])
l = [loss(X, y, b)]
for step in range(n_step):
b -= lr*grad_func(X, y, b)
l.append(loss(X, y, b))
return b, lNow we run the optimization with a variety of different learning rates. For illustration purpose we will only run a few steps.
betas = {}
losses = {}
for lr in [.001, .01, .05, .1, 1]:
betas[lr], losses[lr] = gd_optimize(X, y, grad_func, lr=lr, n_step=15)
for r, b in betas.items():
print("Learning Rate {:5} | Estimate: {}".format(r, b))The result suggests that a learning rate of 1 is too large for our problem. The gradient explodes which make our solution diverge. And a lower learning rate in general converges slower to the optimum. Number of examples used to calculate the gradient also will affect the convergence behavior. In general if the sample size is too small a smaller learning rate should be used to avoid gradient explosion.
This can be more clearly seen if we plot the trace of our training losses:
ofile = "plots/toy_reg_loss_compare.html"
pdata = [go.Scatter(x=np.arange(len(losses[lr])), y=losses[lr], name=lr)
for lr in losses.keys()]
p = plot_offline(
ofile=ofile,
data=pdata,
layout=go.Layout(
title="Trace of Training Loss with Various Learning Rates",
xaxis=dict(title="Step"),
yaxis=dict(title="Loss")))If the loss doesn’t decrease over training iteration, it is a signal that something is wrong with our model.
Unlike our toy implementation, in modern implementation of any numerical optimizer there will be a lot of techniques to do the best to avoid convergence failure. But it is the model developer’s responsibility to diagnose the training behavior before anything is delivered to the stakeholder. Checking the dynamics of loss is usually the first and quick step to examine whether the training task is functioning as expected.
The vanilla gradient descent optimizer we just implemented has one issue. Since for each update it needs to traverse over the entire dataset, it becomes too slow when it comes to large dataset.
Stochastic gradient descent is to overcome this issue. Instead of calculate the gradient using the entire dataset, we can use only one example per update. Each step will be less precise but statistically the final result should be consistent.
Here we introduce the term epoch: One epoch is for the optimizer to traverse the entire dataset once. Number of epochs can be considered as another hyper-parameter of a model.
Here comes the minimum SGD implementation:
def sgd_optimize(X, y, grad_func, lr=.01, n_epoch=10):
b = np.random.normal(size=X.shape[1])
l = [loss(X, y, b)]
out = b.copy()
for epoch in range(n_epoch):
# Shuffle the dataset before each epoch.
sid = np.random.permutation(X.shape[0])
for i in sid:
b -= lr * grad_func(X[None,i,:], y[i], b)
l.append(loss(X[None,i,:], y[i], b))
# Trace the result per epoch.
out = np.row_stack([out, b])
return out, l
sgd_beta, sgd_loss = sgd_optimize(X, y, grad_func, lr=.01, n_epoch=100)
for epoch in [1, 5, 10, 50, 100]:
print("Training Epochs {:4} | Estimate: {}".format(epoch, sgd_beta[epoch]))SGD won’t be as precise as vanilla GD but for large scale application it can reduce considerable amount of computing time (sometimes from infeasible to feasible). In our simple problem indeed just 1 epoch can have a good approximation already.
The trace of (per-instance) loss for the first 500 updates:
ofile = "plots/toy_reg_sgd_loss.html"
p = plot_offline(
ofile=ofile,
data=[go.Scatter(x=np.arange(len(sgd_loss[:500])), y=sgd_loss[:500])],
layout=go.Layout(
title="Trace of SGD Training Loss",
xaxis=dict(title="Step"),
yaxis=dict(title="Loss")))Since the loss is calculated on a per-instance basis, it will fluctuate but with a decreasing trend if nothing went wrong about the optimization process.
In sklearn there are dedicated SGD classes for a variety of learning algorithms. For a linear regression model with SGD solver:
To reduce the noise in SGD we can modify it by replacing 1 training example with a batch of examples in a single gradient update. Here’s such implementation of a batch gradient descent optimizer2:
def gd_batch_optimize(X, y, grad_func, lr=.01, n_epoch=10, batch_size=64):
b = np.random.normal(size=X.shape[1])
l = [loss(X, y, b)]
for epoch in range(n_epoch):
# Shuffle the dataset before each epoch.
sid = np.random.permutation(X.shape[0])
Xs = X[sid,:]
ys = y[sid]
i = 0
n_step = int(np.ceil(X.shape[0] / batch_size))
for step in range(n_step):
Xb = Xs[i:i+batch_size,:]
yb = ys[i:i+batch_size]
b -= lr*grad_func(Xb, yb, b)
l.append(loss(Xb, yb, b))
i += batch_size
return b, l
batch_beta, batch_loss = gd_batch_optimize(X, y, grad_func, n_epoch=100, batch_size=64)
print(batch_beta)For this simple problem batch size doesn’t have any important impact given enough number of training epochs:
for bsize in [8, 16, 32, 64, 128, 256]:
print("Batch Size: {:4} | Estimate: {}".format(
bsize,
gd_batch_optimize(X, y, grad_func, n_epoch=100, batch_size=bsize)[0]))Batch optimizer is currently the best practice of training neural networks in large scale application. The batch size depends on the actual application but usually ranges from 8 (mini) to 1024 (large).
A logistic regression models the outcome probablistically:
\[ \begin{equation} \label{eq:logit} P(Y = 1) = \frac{1}{1 + e^{-X\beta}}, \end{equation} \]
Here the sigmoid function \(s(t) = \frac{1}{1 + e^{-t}}\) is used to transform a real number into probability space \([0, 1]\).
We can interpret the model as a linear model in log-odds. Assuming Y is binary and take a value of 0 or 1, the odds of \(Y = 1\) is defined as \(\frac{P(Y = 1)}{P(Y = 0)} = \frac{P(Y = 1)}{1 - P(Y = 1)}\). We can re-arrange the logistic model equation:
\[ \begin{aligned} \ln \Bigg[ \frac{P(Y = 1)}{1 - P(Y = 1)} \Bigg] &= \ln \Bigg[ \frac{\frac{1}{1 + e^{-X\beta}}}{\frac{e^{-X\beta}}{1 + e^{-X\beta}}} \Bigg] \\ &= \ln(1 + e^{-X\beta}) - \ln e^{-X\beta}(1 + e^{-X\beta}) \\ &= \ln(1 + e^{-X\beta}) - \ln e^{-X\beta} - \ln(1 + e^{-X\beta}) \\ &= - \ln e^{-X\beta} \\ &= X\beta. \end{aligned} \]
That is, the model weights are linear in the log-odds of our target outcome.
When the probability is transformed into odds, the range is transformed from \([0,1]\) to \([0,\infty]\).
# We plot only up to .95 since the value of odds will grow unboundedly.
prob = np.linspace(.001, .95, num=100)
odds = prob / (1 - prob)
log_odds = np.log(odds)
ofile = "plots/prob_odds.html"
p = plot_offline(
ofile=ofile,
data=[go.Scatter(x=prob, y=odds)],
layout=go.Layout(
title="Event Probabilities to Odds",
xaxis=dict(title="Probability"),
yaxis=dict(title="Odds")))If we further transform odds to log-odds, the range is transformed from \([0,\infty]\) to \([-\infty,\infty]\).
ofile = "plots/prob_odds.html"
p = plot_offline(
ofile=ofile,
data=[go.Scatter(x=odds, y=log_odds)],
layout=go.Layout(
title="Event Odds to Log-Odds",
xaxis=dict(title="Odds"),
yaxis=dict(title="Log-Odds")))Put it together is the effect of the sigmoid function:
# Note that instead of coding s = 1 / (1 + np.exp(-t)),
# which is numerically unstable,
# we should use a math trick to make it stable.
def sigmoid(t):
"""Numerically stable sigmoid."""
return np.exp(-np.logaddexp(0, -t))
t = np.linspace(-10, 10, num=100)
ofile = "plots/sigmoid.html"
p = plot_offline(
ofile=ofile,
data=[go.Scatter(x=t, y=sigmoid(t))],
layout=go.Layout(
title="A Sigmoid Function",
xaxis=dict(title="Raw Value"),
yaxis=dict(title="Probability")))How do we solve for the model weights \(\beta\) in equation \(\eqref{eq:logit}\)? In the linear regression model we solve for the weights by minimizing the mean squared error. In logistic regression model we need to define measurement for modeling error as well.
Put it differently, we’d like to calculate the distance between our predicted probability and the real event label distribution (called the empirical distribution). In information theory the cross entropy is used to measure the distance between two probability distribution.
The entropy of a discrete probability distribution is defined as:
\[ H(p) = - \sum_i p_i \log_2p_i, \]
where \(p_i\) is the probability for event \(i\). It measures the uncertainty of a stochastic event.
Take a coin flip event as example. We plot the entropy value at different value of coin bias (the probability of having a head instead of a tail.)
def entropy_binary(p):
return -p*np.log2(p) - (1 - p)*np.log2(1 - p)
prob = np.linspace(.01, .99, num=100)
ofile = "plots/entropy.html"
p = plot_offline(
ofile=ofile,
data=[go.Scatter(x=prob, y=entropy_binary(prob))],
layout=go.Layout(
title="Entropy of a Binary Event (Coin Flip)",
xaxis=dict(title="Probability of a Head"),
yaxis=dict(title="Entropy")))It is obvious that the entropy of this event is maximized when the probability of flipping a head is exactly 0.5. At this level the event has a highest level of uncertainty in a sense that it is the most difficult case to predict the outcome of a flip.
Now move on to cross entropy. Cross entropy between two discrete probability distribution \(p\) and \(q\) over the same support is defined as:
\[ H(p, q) = - \sum_i p_i \log_2q_i. \]
For our logistic regression model, distribution \(p\) is the empirical distribution (label distribution) and distribution \(q\) is our model predicted distribution. Denote \(q_i = P(y_i = 1)\) for the prediction of \(i\)-th example. The cross-entropy-loss of our model hence can be written as:
\[ \mbox{Cross-Entropy Loss} = - \frac{1}{N}\sum_i^N \bigg[ y_i\log_2q_i + (1 - y_i)\log_2(1 - q_i)\bigg], \]
where \(y_i \in \{0,1\}\) is the \(i\)-th binary training label and \(P(y_i = 1)\) the model prediction for the \(i\)-th example out of \(N\) total training examples. This is the mean value of cross entropy for each training example.
Since \(q_i = P(y_i = 1)\) is expressed by our model equation \(\eqref{eq:logit}\), now we can apply gradient descent to the loss function to find out the optimum model weights that minimize the cross-entropy loss.
Before we implement our optimizer for a logistic regression model, we demonstrate that minimizing cross entropy is indeed equivalent to maximizing data likelihood.
The data likelihood is just the product of all predicted probabilities for individual example, provided that they are independent.
\[ \mbox{likelihood} = \prod_i^N q_i^{y_i}(1 - q_i)^{1 - y_i}. \]
The maximum likelihood estimator will try to maximize the log of the likelihood, which is:
\[ \begin{aligned} \mbox{log-lik} &= \log_2 \prod_i^N q_i^{y_i}(1 - q_i)^{1 - y_i} \\ &= \sum_i^N \log_2 q_i^{y_i}(1 - q_i)^{1 - y_i} \\ &= \sum_i^N \bigg[ y_i\log_2 q_i + (1 - y_i)\log_2(1 - q_i) \bigg]. \end{aligned} \]
By taking the average as well, the negative log-likelihood is exactly the cross entropy.
Though we didn’t discuss this linkage in linear regression, if we assume the error term in the model is independently distributed Normally, the OLS solution to the model weights (equation \(\eqref{eq:ols}\)) is indeed also the MLE solution.
Now let’s also implement a batch gradient descent optimizer for a logistic regression model. We will use the same toy example and create additional random binary labels for this exercise.
The gradient function must be derived with respect to model weights. To simplify the math we replace log of base 2 with natural log (with no impact on the optimization), and we take advantage of the fact that the derivative of the sigmoid function is:3
\[ \begin{aligned} \frac{ds(t)}{dt} &= \frac{d\frac{1}{1 + e^{-t}}}{dt} \\ &= -(1 + e^{-t})^{-2} \cdot (- e^{-t}) \\ &= \frac{e^{-t}}{(1 + e^{-t})^2} \\ &= \frac{1}{1 + e^{-t}} \cdot \frac{e^{-t}}{1 + e^{-t}} \\ &= s(t) \cdot [1 - s(t)]. \end{aligned} \]
The gradient in our univariate model w.r.t. the bias term will be:4
\[ \begin{aligned} \frac{\partial\mbox{LogLoss}}{\partial\beta_0} &= - \frac{1}{N}\sum_i \bigg[ y_i\frac{\partial\ln q_i}{\partial\beta_0} + (1 - y_i)\frac{\partial\ln(1 - q_i)}{\partial\beta_0}\bigg] \\ &= - \frac{1}{N}\sum_i \bigg[ y_i\frac{1}{q_i}\frac{\partial q_i}{\partial\beta_0} + (1 - y_i)\frac{1}{1 - q_i}\frac{\partial(1 - q_i)}{\partial\beta_0}\bigg] \\ &= - \frac{1}{N}\sum_i \bigg[ y_i\frac{1}{q_i}\frac{\partial (\beta_0 + \beta_1x_i)}{\partial\beta_0}\frac{\partial q_i(t)}{\partial t} - (1 - y_i)\frac{1}{1 - q_i}\frac{\partial (\beta_0 + \beta_1x_i)}{\partial\beta_0}\frac{\partial q_i(t)}{\partial t}\bigg] \\ &= - \frac{1}{N}\sum_i \bigg[ y_i\frac{1}{q_i}q_i(1 - q_i) - (1 - y_i)\frac{1}{1 - q_i}q_i(1 - q_i)\bigg] \\ &= - \frac{1}{N}\sum_i \bigg[ y_i(1 - q_i) - (1 - y_i)q_i \bigg] \\ &= - \frac{1}{N}\sum_i (y_i - q_i). \end{aligned} \]
Similary for the weight:
\[ \frac{\partial\mbox{LogLoss}}{\partial\beta_1} = - \frac{1}{N}\sum_i (y_i - q_i)x_i. \]
Now the Python code for batch gradient descent with log-loss:
def logloss(X, y, b):
logloss_pos = y * np.log(sigmoid(X.dot(b)))
logloss_neg = (1 - y) * np.log(1 - sigmoid(X.dot(b)))
return - (logloss_pos + logloss_neg).mean()
def logit_grad_func(X, y, b):
return - (y - sigmoid(X.dot(b))).dot(X) / X.shape[0]
def cross_entropy_gd_optimize(X, y, lr=.01, n_epoch=10, batch_size=64):
b = np.random.normal(size=X.shape[1])
l = [logloss(X, y, b)]
for epoch in range(n_epoch):
# Shuffle the dataset before each epoch.
sid = np.random.permutation(X.shape[0])
Xs = X[sid,:]
ys = y[sid]
i = 0
n_step = int(np.ceil(X.shape[0] / batch_size))
for step in range(n_step):
Xb = Xs[i:i+batch_size,:]
yb = ys[i:i+batch_size]
b -= lr*logit_grad_func(Xb, yb, b)
l.append(logloss(Xb, yb, b))
i += batch_size
return b, l
log_beta, log_loss = cross_entropy_gd_optimize(X, y2, lr=.01, n_epoch=100)
print(log_beta)ofile = "plots/toy_log_loss.html"
pdata = [go.Scatter(x=np.arange(len(log_loss)), y=log_loss)]
p = plot_offline(
ofile=ofile,
data=pdata,
layout=go.Layout(
title="Trace of Training Log-Loss",
xaxis=dict(title="Step"),
yaxis=dict(title="Loss")))Comparing to linear regression, a logistic regression model is harder to converge. Since our naive implementation does not do convergence diagnostics, let’s use R’s built-in glm function which use Newton’s method (a 2nd-order optimizer utilizing not only 1st-order but also 2nd-order derivatives) to check the estimation result of our toy example:
Now increase both the learning rate and training epochs of our naive gradient descent optimizer:
Seems better.
Or we can use the Python package statsmodels (which also uses a 2nd-order optimizer by default) to check our result:
Or the sklearn result:
from sklearn.linear_model import LogisticRegression
# We need to set param C to an arbitrarily large number to suppress regularization
# since our naive approach doesn't implement any regularization.
# We set fit_intercept=False since our design matrix already contains intercept.
print(LogisticRegression(C=1e16, fit_intercept=False, solver="liblinear").fit(X, y2).coef_)In the previous section we implement a simple gradient descent optimizer by manually derive the functional form of gradient on our own. This could be troublesome if our model becomes more and more complicated, as in the case of a deep neural net.
Automatci differentiation is a programming technique to calculate the gradient of any given function. One of the most popular library for this purpose is TensorFlow.
Let’s use tensorflow to implement our simple gradient descent optimizer again. But this time we will NOT explicitly derive the gradient function. Instead, we will only specify the target function which is just the loss function of our model.
import tensorflow as tf
tf.enable_eager_execution()
# Use tensor to represent our data.
# Note that we need to be very specific about dtype/shape of our tensors.
X_tf = tf.convert_to_tensor(X, dtype=tf.float32)
beta_tf = tf.reshape(tf.convert_to_tensor(beta, dtype=tf.float32), (2,1))
e_tf = tf.reshape(tf.convert_to_tensor(e, dtype=tf.float32), (N, 1))
y_tf = tf.matmul(X_tf, beta_tf) + e_tf
y2_tf = tf.reshape(tf.convert_to_tensor(y2, dtype=tf.float32), (N, 1))
def loss_mse_tf(X, y, beta):
y_hat = tf.matmul(X, beta)
return tf.reduce_mean(input_tensor = (y_hat - y)**2)
def logloss_tf_bad(X, y, beta):
# This can suffer from numerical stability issue.
logloss_pos = y * tf.log(tf.sigmoid(tf.matmul(X, beta)))
logloss_neg = (1 - y) * tf.log(1 - tf.sigmoid(tf.matmul(X, beta)))
return - tf.reduce_mean(logloss_pos + logloss_neg)
def logloss_tf(X, y, beta):
# tf.sigmoid is not numerically stable.
# But there is a convenient function to do the cross entropy calculation stably.
s = tf.nn.sigmoid_cross_entropy_with_logits(labels=y, logits=tf.matmul(X, beta))
return tf.reduce_mean(s)
def gd_optimize_tf(X, y, loss_func, lr=.01, n_step=100):
grad_func_tf = tf.contrib.eager.gradients_function(loss_func, params=[2])
beta = tf.random_normal((2, 1))
for step in range(n_step):
grad = grad_func_tf(X, y, beta)[0]
beta -= lr * grad
return beta.numpy()
# For the MSE loss problem.
print(gd_optimize_tf(X_tf, y_tf, loss_func=loss_mse_tf, n_step=3000))
# For the cross entropy loss problem.
print(gd_optimize_tf(X_tf, y2_tf, loss_func=logloss_tf, lr=.5, n_step=5000))In the above coding example one should realize that we no longer need to hardcode the functional form of the gradients. Instead we just plug-in the loss function and let tensorflow to do the gradient calculation for us.
Of course in actual development we will use higher-level APIs to implement our model, where the entire optimization process is abstracted away from the application code.
Both linear regression and logistic regression can be considered as simple additive model of the form:
\[ \hat{y} = \Phi\bigg(\sum_{i=1}^Pw_ix_i\bigg), \]
where \(P\) is the number of features used, \(x_i\) is the \(i\)-th feature, and \(\Phi(\cdot)\) is a function applied to the output. In linear regression \(\Phi(\cdot)\) is simply an identity function. In logistic regression \(\Phi(\cdot)\) is the standard sigmoid function.
Now consider there is a way to ensemble multiple such additive models together to generate a potentially better and more sophisticated model. We use the following diagram for illustration.
In the above diagram, \(Y_{11}\) is a single additive model
\[ Y_{11} = \Phi(W_{111}X_1 + W_{211}X_2 + W_{311}X_3). \]
Similarly \(Y_{12}\) is another such model (with the same input feature set but different model weights)
\[ Y_{12} = \Phi(W_{112}X_1 + W_{212}X_2 + W_{312}X_3). \]
Model \(Y_2\) is yet another additive model but takes the output of the above two models:
\[ Y_2 = \Phi(W_{1121}Y_{11} + W_{1221}Y_{12}). \]
The above setup is a simple architecture of a neural network model, with only one hidden layer of two neurons. (We also ignore the constant/bias term in each layer for simplicity.) A neuron is simply an additive model with a so-called activation function \(\Phi(\cdot)\) to transform the output from any real number into a scaled signal.
One now can easily realize that a logistic regression model could be viewed as a degenerated neural network model with single neuron and using sigmoid as the activation function. And a linear regression model is a degenerated neural network model with single neuron and without an activation function.
Why do we need the activation function? In the above neural network model in the absence of activation function the final output model \(Y_2\) will degenerate into a simple linear model. We can use simpler notations to demonstrate this:
\[ \begin{aligned} y_1 &= ax + b, \\ y_2 &= cx + d, \\ y_3 &= e + fy_1 + gy_2 \\ &= e + f(ax + b) + g(cx + d) \\ &= \underbrace{(e + fb + gd)}_\text{Bias} + \underbrace{(fa + gc)}_\text{Weight}x. \end{aligned} \]
Without activation function, no matter how many neurons or layers we design for our model, it eventually reduces to a simple linear model. With the activation function applied to each neuron, the model becomes non-linear and hence can handle much more complicated patterns hidden behind the data.
Some popular activation functions:
To solve for model weights in a neural network, we use a technique called backpropagation which is essentially an iterative process of gradient descent.
To simplify notation we assume each neuron is simply a univariate model. Consider the following minimum architecture:
Mathematically:
\[ \begin{aligned} \hat{y} &= \Phi(b_1 + w_1z) \\ &= \Phi(b_1 + w_1\Phi(b_0 + w_0x)), \end{aligned} \]
where
\[ z = \Phi(b_0 + w_0x), \]
and
\[ \Phi(t) = \frac{1}{1 + e^{-t}}. \]
Though it may not be very meaningful to use MSE loss when the output layer is applied with an activation function, we can still do it for educational purpose.
\[ \mbox{MSE-Loss} = \frac{1}{N}\sum_i^N (y_i - \hat{y}_i)^2. \]
Firstly we take the derivative w.r.t. the weight in the last layer:
\[ \begin{aligned} \frac{\partial \mbox{MSE-Loss}}{\partial w_1} &= - \frac{1}{N}\sum_i (y_i - \hat{y}_i)\frac{\partial \hat{y}_i}{\partial w_1} \\ &= - \frac{1}{N}\sum_i (y_i - \hat{y}_i) \underbrace{ \frac{\partial t}{\partial w_1}\frac{\partial \Phi(t)}{\partial t} }_{t = b_1 + w_1\Phi(b_0 + w_0x)}\\ &= - \frac{1}{N}\sum_i (y_i - \hat{y}_i) \cdot \Phi(t)(1 - \Phi(t)) \cdot z_i. \end{aligned} \]
Similarly for the bias in the last layer:
\[ \frac{\partial \mbox{MSE-Loss}}{\partial b_1} = - \frac{1}{N}\sum_i (y_i - \hat{y}_i) \cdot \Phi(t)(1 - \Phi(t)). \]
Now move on to the bias and weight in the first layer:
\[ \begin{aligned} \frac{\partial \mbox{MSE-Loss}}{\partial w_0} &= - \frac{1}{N}\sum_i (y_i - \hat{y}_i)\frac{\partial \hat{y}_i}{\partial w_0} \\ &= - \frac{1}{N}\sum_i (y_i - \hat{y}_i) \underbrace{ \frac{\partial t}{\partial w_0}\frac{\partial \Phi(t)}{\partial t} }_{t = b_1 + w_1\Phi(b_0 + w_0x)}\\ &= - \frac{1}{N}\sum_i (y_i - \hat{y}_i) \cdot \Phi(t)(1 - \Phi(t)) \cdot \underbrace{ \Phi(k)(1 - \Phi(k)) }_{k = b_0 +w_0x} \cdot w_1 \cdot x_i, \\ \frac{\partial \mbox{MSE-Loss}}{\partial b_0} &= - \frac{1}{N}\sum_i (y_i - \hat{y}_i) \cdot \Phi(t)(1 - \Phi(t)) \cdot \Phi(k)(1 - \Phi(k)). \end{aligned} \]
One can clearly see there is a linkage between the derivative of the weights in consecutive layers:
\[ \begin{aligned} \frac{\partial \mbox{MSE-Loss}}{\partial w_1} &= - \frac{1}{N}\sum_i (y_i - \hat{y}_i) \cdot \Phi(t)(1 - \Phi(t)) \cdot z_i, \\ \frac{\partial \mbox{MSE-Loss}}{\partial w_0} &= - \frac{1}{N}\sum_i (y_i - \hat{y}_i) \cdot \underbrace{\Phi(t)(1 - \Phi(t))}_{\hat{y}_i(1 - \hat{y}_i)} \cdot \underbrace{\Phi(k)(1 - \Phi(k)) \cdot w_1 \cdot x_i}_{\frac{\partial w_1z}{\partial w_0} = \frac{\partial w_1\Phi(b0 + w_0x)}{\partial w_0}}. \end{aligned} \]
Similarly we can derive the gradients for cross-entropy loss.
\[ \begin{aligned} \mbox{LogLoss} &= - \frac{1}{N}\sum_i^N \bigg[ y_i\ln\hat{y}_i + (1 - y_i)\ln(1 - \hat{y}_i)\bigg], \\ \frac{\partial \mbox{LogLoss}}{\partial w_1} &= - \frac{1}{N} \sum_i (y_i - \hat{y}_i)z_i, \\ \frac{\partial \mbox{LogLoss}}{\partial w_0} &= - \frac{1}{N} \sum_i (y_i - \hat{y}_i)\Phi(k)(1 - \Phi(k))w_1x_i.\\ \end{aligned} \]
(Skip the bias terms to save space.)
Now let’s implement the simple neural network model in Python:
def sigmoid(x):
return np.exp(-np.logaddexp(0, -x))
def sigmoid_grad(z):
return z*(1 - z)
class NeuralNetwork:
"""Simple neural network with 1 hidden layer of 4 neurons."""
def __init__(self, x, y):
self.input = x
self.w1 = np.random.rand(self.input.shape[1], 4)
self.w2 = np.random.rand(4, 1)
self.y = y
self.output = np.zeros(y.shape)
def feedforward(self):
self.layer1 = sigmoid(np.dot(self.input, self.w1))
self.output = sigmoid(np.dot(self.layer1, self.w2))
def backprop(self, lr):
raise NotImplementedError
def train(self, lr, n_step):
for step in range(n_step):
self.feedforward()
self.backprop(lr=lr)
yhat = ",".join([str(y) for y in np.squeeze(self.output)])
return yhat
class NerualNetworkWithMSELoss(NeuralNetwork):
def backprop(self, lr):
# We ignore all constant terms in the gradient.
self.dy = self.y - self.output
dw2 = - np.dot(self.layer1.T, (self.dy*sigmoid_grad(self.output)))
dw1 = - np.dot(self.input.T,
(np.dot(self.dy*sigmoid_grad(self.output),
self.w2.T)*sigmoid_grad(self.layer1)))
self.w1 -= lr * dw1
self.w2 -= lr * dw2
class NerualNetworkWithLogLoss(NeuralNetwork):
def backprop(self, lr):
# We ignore all constant terms in the gradient.
self.dy = self.y - self.output
dw2 = - np.dot(self.layer1.T, self.dy)
dw1 = - np.dot(self.input.T,
np.dot(self.dy, self.w2.T)*sigmoid_grad(self.layer1))
self.w1 -= lr * dw1
self.w2 -= lr * dw2Here we will use the hello-world example of artificial neural network: The XOR problem. The XOR logical outcome, besides extremely simple, is not linearly separable. So it serves as a good example of showcasing neural networks’ non-linearity.
The input data is simply combination of two binary switches. (For completeness we also include a constant term which always evaluate to 1 as the first feature.) The output data is the XOR result.
Xnn = np.array(
[[1, 0, 0],
[1, 1, 0],
[1, 0, 1],
[1, 1, 1]])
ynn = np.array([[0], [1], [1], [0]])
print(Xnn) # Input features: Bias + two binary switches.
print(ynn) # XOR output.Now let’s see if our simple neural network model can learn the XOR pattern.
nn_mse = NerualNetworkWithMSELoss(Xnn, ynn)
nn_mse.train(lr=1, n_step=3000) # The predicted XOR outcome.
print(nn_mse.w1) # Estimated 1st-layer model weights
print(nn_mse.w2) # Estimated 2nd-layer model weightsFor this simple example, both MSE and cross-entropy loss can work fine to figure out the XOR pattern. (Strictly speaking cross-entropy loss performs better.)
Regularization is a technique to mitigate overfitting. It is not particular to neural networks but general to all machine learning models.
A regularization is a constraint added onto the original optimization problem (loss minimization). Consider model weights as a real vector, we usually use the norm of the vector to constraint its size in the optimization. The norm is a measure of size of a vector. The general idea is to add a penalty to the target function (which is a loss function) such that it discourages using large weights (which increase the size of vector) to achieve the minimization goal.
In general, a MSE loss minimization with a p-norm regularization for a linear model can be written (in matrix form) as:
\[ \begin{aligned} \hat{y} &= X\mathrm{B}, \\ \operatorname*{argmin}_\beta \mbox{Loss} &= \frac{1}{N}\bigg[ \underbrace{(y - X\mathrm{B})^T(y - X\mathrm{B})}_\text{sum of squared errors} + \underbrace{\lambda\Vert\mathrm{B}\Vert_{p}}_\text{p-norm regularization} \bigg], \end{aligned} \]
with a model weight vector \(\mathrm{B} = [\beta_1, \beta2, ..., \beta_k]\). The parameter \(\lambda\) is a new hyper-parameter introduced by the regularization.
Or put it in a scalar form:
\[ \operatorname*{argmin}_\beta \mbox{Loss} = \frac{1}{N}\sum_{i=1}^N\bigg( y_i - \mathrm{B}x_i\bigg)^2 + \frac{\lambda}{N}\bigg(\sum_{j=1}^k\vert\beta_j\vert^p\bigg)^{1/p}. \]
Popular choices of \(p\) is \(p = 1\) (L1-Norm) and \(p = 2\) (L2-Norm). Notably, \(p\) can be any non-integer real value and would have even better result than integer-norm. But their computational difficulty make them much less desirable in practice.
In practice, we also ignore the root in the norm for computational simplicity. So the problem indeed becomes
\[ \operatorname*{argmin}_\beta \mbox{Loss} = \frac{1}{N}\sum_{i=1}^N\bigg( y_i - \mathrm{B}x_i\bigg)^2 + \frac{\lambda}{N}\bigg(\sum_{j=1}^k\vert\beta_j\vert^p\bigg). \]
There is also a so-called elasticnet regularization which combines both L1 and L2 norms together to give potentially better result. However it requires a lot more computing resources and hence become much less feasible in a neural network model.
The L2-Norm shrinks the size of model weights. Those who contribute less to loss minimization will shrink more.
Back to our univariate linear model \(\hat{y}_i = \beta_0 + \beta_1x_i\) for symolic simplicity. Our MSE loss is
\[ \mbox{MSE-Loss} = \frac{1}{N}\sum_i\big[y_i-(\beta_0 + \beta_1x_i)\big]^2 + \frac{\lambda}{N}(\beta_0^2 + \beta_1^2). \]
Since the penaty is additive, our gradient solution is extremely easy:
\[ \begin{aligned} \frac{\partial\mbox{MSE-Loss}}{\partial\beta_0} &= - \frac{2}{N}\sum_i \big[ y_i - (\beta_0 + \beta_1x_i) \big] + \frac{2}{N}\lambda\beta_0, \\ \frac{\partial\mbox{MSE-Loss}}{\partial\beta_1} &= - \frac{2}{N}\sum_i \big[ y_i - (\beta_0 + \beta_1x_i) \big]x_i + \frac{2}{N}\lambda\beta_1. \end{aligned} \]
A linear regression model with l2 regularization is referred to as the ridge regression. Be aware that in implementing the batch optimizer we need to scale down the update from l2 norm by the size of batch in order to balance the size of l2 norm and the original loss.
def mse_grad_l2(X, y, beta):
a = 1 # Hardcode the regularization coefficient for simplicity.
return -2 * ((y - X.dot(beta)).dot(X)) / X.shape[0] + 2 * a * beta / X.shape[0]
batch_beta, batch_loss = gd_batch_optimize(X, y, mse_grad_l2, n_epoch=100, batch_size=64)
print(batch_beta)The resulting estimates are smaller than the vanilla gradient descent.
One thing worth noting is that when batch gradient descent is used along with regularization, a side effect exists due to the scaling of the batch size on the regularization term.
Experiments with different batch sizes with l2 regularization:
for bsize in [8, 16, 32, 64, 128, 256, 512, 1000]:
print("Batch Size: {:4} | Estimate: {}".format(
bsize,
gd_batch_optimize(X, y, mse_grad_l2, n_epoch=500, batch_size=bsize)[0]))Let’s check the result using sklearn’s Ridge regressor (which doesn’t use gradient descent as its solver):
# We set fit_intercept=False since our design matrix already contains intercept.
# According to the doc the constant 2 in the derivative is not included so we set alpha = 2.
from sklearn.linear_model import Ridge
print(Ridge(alpha=2, fit_intercept=False).fit(X, y).coef_)The result is closed to gradient descent with nearly a full-batch update (i.e., use the entire training data for one update).
Exactly the same logic can apply to logistic regression and is not discussed here to avoid redundancy.
The L1-Norm prefers sparsity in the model weights. That is, it may zero out weights that are not contributing to the loss. While the L2-Norm only makes them arbitrarily small.
MSE loss with l1 norm:
\[ \mbox{MSE-Loss} = \frac{1}{N}\sum_i\big[y_i-(\beta_0 + \beta_1x_i)\big]^2 + \lambda(\vert\beta_0\vert + \vert\beta_1\vert), \]
with gradients:5
\[ \begin{aligned} \frac{\partial\mbox{MSE-Loss}}{\partial\beta_0} &= - \frac{2}{N}\sum_i \big[ y_i - (\beta_0 + \beta_1x_i) \big] + \lambda\frac{\beta_0}{\vert\beta_0\vert}, \\ \frac{\partial\mbox{MSE-Loss}}{\partial\beta_1} &= - \frac{2}{N}\sum_i \big[ y_i - (\beta_0 + \beta_1x_i) \big]x_i + \lambda\frac{\beta_1}{\vert\beta_1\vert}. \end{aligned} \]
A linear regression model with l1 regularization is referred to as the lasso regression. Let’s first check the result using sklearn’s Lasso regressor (which uses a coordinate descent solver):
# We set fit_intercept=False since our design matrix already contains intercept.
from sklearn.linear_model import Lasso
print(Lasso(alpha=1, fit_intercept=False).fit(X, y).coef_)To derive roughly the same result, we need to calibrate the coefficient on regularization in our gradient descent solver:
def mse_grad_l1(X, y, beta):
# a is set to 2 since in sklearn convention the entire loss is scaled by a factor of 1/2
# to cancel out the derivative constant on the squared error term.
a = 2
return -2 * ((y - X.dot(beta)).dot(X)) / X.shape[0] + a * np.sign(beta)
# We use a full-batch gradient descent here.
batch_beta, batch_loss = gd_batch_optimize(X, y, mse_grad_l1, n_epoch=1000, batch_size=1000)
print(batch_beta)We can also check the consistency by calling sklearn’s SGDRegressor:
from sklearn.linear_model import SGDRegressor
print(SGDRegressor(penalty="l1", alpha=1, fit_intercept=False).fit(X, y).coef_)Now let’s purposely add one random feature as a noise into the design matrix and see how the regularization helps to shrink the weight of the noise.
# Introduce a totally irrelevant feature into our linear model.
noise = np.random.uniform(size=N).reshape(N, 1)
Xnoise = np.concatenate([X, noise], axis=1)
print(Ridge(alpha=2, fit_intercept=False).fit(Xnoise, y).coef_)
print(Lasso(alpha=1, fit_intercept=False).fit(Xnoise, y).coef_)As one can see, Lasso regressor is able to completely wipe out the weight on the redundant noise, while Ridge can only shrink it toward zero.
To arrive roughly at the same result of Ridge using our simple gradient descent solver:
And for the Lasso case:
One thing to note: It is in general difficult to arrive at sparse solution (0 weight) for a naive gradient descent implementation such as ours. Some further mathematical tricks like the proximal gradient metohd must be applied in order to produce sparse result under finite iterations.
Finally, again, exactly the same logic can apply to logistic regression and is not discussed further to avoid redundancy.
If we took a Bayesian approach on the weight estimation, then a Normal prior on the weight will result in a maximum a posteriori estimator to have a target function with the l2 regularization term, and a Laplacean prior will result in the l1 term.
from scipy.stats import laplace, norm
xticks = np.linspace(-5, 5, num=1000)
ofile = "plots/reg_prior.html"
p = plot_offline(
ofile=ofile,
data=[go.Scatter(x=xticks, y=laplace.pdf(xticks), name="Laplacean Prior (L1)"),
go.Scatter(x=xticks, y=norm.pdf(xticks), name="Normal Prior (L2)")],
layout=go.Layout(
title="Bayesian Prior on Model Weight",
xaxis=dict(title="Model Weight"),
yaxis=dict(title="Density")))Compared to a Normal prior, a Laplacean put much more density at the center. When both distribution is standardized (with location at 0 and scale at 1) as just plotted above, this can be interpreted as that a Laplacean prior has a stronger preference for a sparse solution (weight = 0).
Another maybe easier way of regularization is the dropout mechanism. It is first introduced in a neural network model but can apply to some other machine learning algorithms as well (such as gradient boosting trees).
One can try change the sample size arbitrarily to see the behavior of the estimator.↩
Sometimes “mini-batch gradient descent” is used to describe SGD with batch size > 1. Here we ignore the redundant wording and just call it batch gradient descent.↩
Here we use the common derivatives: \(\frac{df(x)^n}{dx} = nf(x)^{n-1}f'(x)\) and \(\frac{de^{f(x)}}{dx} = f'(x)e^{f(x)}.\)↩
Here we use the common derivative \(\frac{d\ln f(x)}{dx} = \frac{f'(x)}{f(x)}\) and the chain rule to handle the term \(\frac{\partial\ln q_i}{\partial\beta_0}\) where \(q_i\) is the sigmoid computed for the \(i\)-th example.↩
Note that the derivative \(\frac{du}{dx} = \frac{u}{\vert u \vert}\frac{du}{dx}\) at \(u = 0\) is undefined.↩